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Abstract 

Monte Carlo simulations of systems with a complex action are known to be extremely 
difficult. A new approach to this problem based on a factorization property of distribu- 
tion functions of observables has been proposed recently. The method can be applied 
to any system with a complex action, and it eliminates the so-called overlap problem 
completely. We test the new approach in a Random Matrix Theory for finite density 
QCD, where we are able to reproduce the exact results for the quark number density. 
The achieved system size is large enough to extract the thermodynamic limit. Our 
results provide a clear understanding of how the expected first order phase transition 
is induced by the imaginary part of the action. 
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1 Introduction 



The (Euclidean) action of many interesting systems in fields ranging from condensed matter 
physics to high-energy physics has an imaginary part. Some examples in high-energy physics 
are QCD at finite baryon density, Chern-Simons theories, systems with topological terms (like 
the #-term in QCD) and systems with chiral fermions. While this is not a conceptual problem 
per se, it severely limits the application of Monte Carlo methods, which otherwise might 
provide a powerful tool to understand the properties of these systems from first principles. 
So far there is no general solution to this 'complex action problem'. 

In Ref. a new Monte Carlo approach to systems with a complex action was proposed. 
This method utilizes a simple factorization property of distribution functions of observables. 
Since the property holds quite generally, the approach can be applied to any system with a 
complex action. Most notably, the method eliminates the so-called overlap problem, which 
occurs when one applies the standard re-weighting technique to include the effect of the 
imaginary part. Ultimately we hope that this method will enable us, among other things, 
to explore the phase diagram of QCD at finite baryon density, where interesting phases such 
as a superconducting phase have been conjectured to appear ||. 

As a first step toward achieving this goal we test the new approach in a Random Matrix 
Theory for finite density QCD ||. Random Matrix Theory was originally introduced to 
describe the spectrum of the Dirac operator at zero chemical potential and has been 
studied intensively in the literature. (See Ref. @ for a review). The particular extended 
model we study can be regarded as a schematic model for QCD at finite baryon density. As 
one increases the 'chemical potential', the model undergoes a first order phase transition, 
where the imaginary part of the action plays a crucial role. Since it is solvable even for finite 
matrix size iV ||, it serves as a useful testing ground for simulation techniques for QCD at 
finite density. For instance, the problem with quenched simulations || has been clarified in 
Ref. ||. It was also used to test the so-called Glasgow method [fTD| , arid the source of the 
problems was identified |jTT] . 

In this article we apply the new method to both phases of this model (below and above 
the critical point) and obtain the expectation value of the 'quark number density'. The 
results nicely reproduce the exact results known for finite N. The values of iV that are 
accessible by the new method turn out to be large enough to extract the large- N limit. 
Moreover, our results provide a clear understanding of how the first order phase transition 
is induced by the imaginary part of the action. 

The method [|12| proposed for simulating ^-vacuum like systems can be regarded as a 
special case of the factorization method. A simplified version of the method was sufficient 
because the observable was identical to the imaginary part of the action. The essence of 
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the factorization method is that it avoids the overlap problem by the use of constrained 
simulations. Results for 2d CP 3 and other models are very promising. 

The remainder of this article is organized as follows. In Section [| we define a Random 
Matrix Theory for QCD at finite baryon density and review the known exact results. In 
Section [| we explain the complex action problem associated with the standard re-weighting 
technique. The application of the factorization method to a Monte Carlo study of the 
Random Matrix Model is discussed in Section |j. Section |5] shows the results which nicely 
reproduce the exact values for the quark number density. Section § is devoted to summary 
and discussions. 



2 Random Matrix Theory for finite density QCD 

The Random Matrix Model we study in this article is defined by the partition function 

Z= f dWe~ Nt < w ' w ^ detD , (2.1} 



where W is a N x N complex matrix, and the D is a 2N x 2N matrix given by 

m iW + n 
i + n m 



The parameters m and \x correspond to the 'quark mass' and the 'chemical potential', re- 
spectively. The size of the matrix can be thought of as the total number of 'low- lying' modes 
for given total volume. Since the density of these modes is taken to be unity, N can be 
interpreted as the volume of space time. This model has the global symmetries of the Dirac 
operator of QCD at nonzero baryon chemical potential, where the matrix W in ( |2.2| ) is a 
complicated function of the background gauge field. Thus the above model can be thought 
of as a schematic model for QCD at finite baryon density, where the path integral over the 
gauge field is simply replaced by the Gaussian integral over W. Interesting observables are 
the 'chiral condensate' and the 'quark number density' defined by 

2 = ^'■'(■D- 1 ) (2.3) 

The model was first solved in the large- N limit H. Later it was noticed that the model 
can be solved even for finite iV ||. Throughout this paper, we consider the massless case 
(m = 0) for simplicity. Then the partition function can be expressed as 
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1+ m 7(^ + 1, ft) 



(2.5) 
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Figure 1: The exact result ( |2.8| ) for the 'quark number density' (u) is plotted as a function 
of the 'chemical potential' /i for N = 8, 16, 32, 64, 128. In the N — > oo limit, the function 
develops a discontinuity at \i = /i c = 0.527 



where k = —Nfi 2 and 7(71,2;) is the incomplete 7-function defined by 



7(n,x) 



(2.6) 



From this one obtains the vacuum expectation value (VEV) of the quark number density as 

1 d 



2N dfi 



\nZ(fi) 



1 + 



N - 

k e 



(2.7) 
(2.8) 



(-l) iV+1 iV+7(A^ + 1,k)_ 

In Fig. |I] we plot (v) as a function of the chemical potential \x for = 8, 16, 32, 64, 128. The 
large-iV limit of this formula can be taken easily by applying the saddle-point analysis to 
the incomplete 7-function. We obtain 



lim (v) 

N^QO 



—fi for /i < fi c 
1/ ji for \i > fi c , 



(2.9) 



where fi c is the solution to the equation 1 + /i 2 + ln(/i 2 ) = 0, and its numerical value is given 
by \i c = 0.527 • • •. We find that the quark number density (u) has a discontinuity at \x = /x c . 
Thus the schematic model reproduces qualitatively the first order phase transition expected 
to occur in 'real' QCD at nonzero baryon density 
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3 The complex action problem 



In this section we describe the complex action problem that appears in standard Monte 
Carlo studies of the model (|2.1|). Let us first rewrite (|2.1|) as 

dWe~ So+ir , (3.1) 

where we have introduced So and T by 

S = Nti:(W^W) -\n\detD\ (3.2) 
det£> = e ir | det D| . (3.3) 

In this form it becomes manifest that the system has a complex action, where the problematic 
imaginary part T is given by the phase of the fermion determinant. Since the weight e 0+ * r 
in (|3.1|) is not positive definite, we cannot regard it as a probability density. Hence it seems 
difficult to apply the idea of standard Monte Carlo simulations, which reduces the problem of 
obtaining VEVs to that of taking an average over an ensemble generated by the probability 
density. One way to proceed is to apply the reweighting method and rewrite the VEV (v) 
as 

w = ^ • (3-4) 

where the symbol ( ■ )o denotes a VEV with respect to the so-called phase quenched partition 
function 

Z = J dWe- Ntl( - w ^ | det£>| = J dW e~ So . (3.5) 



Since the system (|3.5| ) has a positive definite weight, the VEV ( • )o can be evaluated by 
standard Monte Carlo simulations. However, the fluctuations of the phase V in (|3.4j) grows 



linearly with the size of the matrix D, which is of O(N). Due to huge cancellations, both the 
denominator and the numerator of the r.h.s. of (|3.4j ) vanish as e~ const iV as N increases, while 
the 'observables' e* r and ve lV are of 0(1) for each configuration. As a result, the number 
of configurations required to obtain the VEVs with some fixed accuracy grows as e constJV 
(We remind the reader that N can be considered as the volume of space time). This is 
the notorious 'complex action problem' (or rather the 'sign problem', as we see below). See 
Refs. ]T3| for simulation results for 'real' finite density QCD obtained by this reweighting 
technique. 

In fact we may simplify the expression (|3.4j) slightly by using a symmetry. We note 
that the fermion determinant det D, as well as the observable v, becomes complex conjugate 
under the transformation 

W h-> —W , (3.6) 
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while the Gaussian action remains invariant. From this we find that 



<") = (m)+i (vi) (3.7) 
(^psT), M = i <«*gh t (3 . 8) 

(cosr) (cosr) 

where z/r and z/j denote the real part and the imaginary part of u, respectively. In Eq. 
( |3.8|) the problem takes the form of the 'sign problem', since cosT and sinT flip their sign 



violently as a function of the configuration W. Note that both terms in the r.h.s. of (|3.7| ) 
are real, meaning in particular that their sum (u) is also real. 

The model (|3.5|) is solvable in the large- N limit ||. For m = one obtains 

lim M, = ( T^] (3-9) 

tv^oo ' [ l/A* for > 1 . v 7 

In this case the VEV of the quark number density is a continuous function of the chemical 
potential /i unlike in (|2.9|) . This remarkable difference between ( [2.9D and (|3.9| ) is precisely 
due to the imaginary part Y of the action. The two results agree trivially at fi — 0, but 
interestingly they also agree at \i > 1. This is because the eigenvalues of W are located 



inside the complex unit circle [11 . If // > 1, the effect of the fermion determinant is 1/iV 
suppressed, and the quenched approximation, as well as the phase quenched approximation, 
becomes exact in the thermodynamic limit. Note also that the symmetry under (|3.6|) implies 



fa)o = ; Wo = Mo. (3.10) 

We can also predict (i/r) and (z/[) separately for the unquenched model in the large N 
limit. For that we note that v{— /i) = —v(pL)*. Therefore we have 

W = \{("<J*)) ~ M-A*)>}, (^i) = ^{Mm)) + (K-A*))} , (3.H) 

where in the calculation of the sign of the chemical potential in the fermion deter- 

minant is not changed. In calculating (z/(/i)), the fermion determinant det D(fi) cancels the 
poles of This is responsible for the difference between quenched and unquenched re- 

sults for fi < 1. However, the cancellation of the poles does not occur in calculating {v{—fi)). 
Therefore we expect that (u(—fi)) coincides with the corresponding quenched result (v(— //)) 
even for /i < 1. This implies in particular that for \x < /i c we obtain the results (z/r) = and 
i (yi) = —fi in sharp contrast to the quenched result (|3.10| ). For fj, > 1, on the other hand, 
both (vr) and (z/j) agree with the corresponding quenched results. These results are indeed 
observed in our Monte Carlo simulations, which shall be discussed in Section |5]. 
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4 The factorization method 
4.1 The basic formulae 

In this section, we explain how the factorization method can be used to obtain the VEVs 
(z/ R ) and (ui). The fundamental objects of the method are the distribution functions 

p R (x) d = f (S(x-Re(u))) and p l (y) d = f (5(y - Irm»)) (4.1) 

associated with the complex valued v . In a unified notation that will be used below these 
equations can be rewritten as 

Pi {x) = (6(x - Ui)) i = R,I. (4.2) 
In terms of these functions, the VEVs of V{ can be expressed as 

/oo 

dxxpi(x) . (4.3) 
-oo 

What is essential for the method are the constrained partition functions 

Zi{x) = J dWe~ So 5{x - u { ) (4.4) 

and the average of the phase e* r with respect to these partition functions 

<Pi(x) = (e ir ) M . (4.5) 
If we also introduce the distribution functions for the phase quenched model by 

p 4 {0) (x) d = f (*(s-i/ 4 )> , (4.6) 
the distribution functions for the unquenched model trivially factorize as 

p l (x) = ^pf\x) Vl (x) i = R, I, (4.7) 
where the normalization constant C is given by 

C^<e* r >o. (4.8) 

Plugging fl4.7| ) into (f4.3|), we get 

1 f°° 

(»i) = cj dxxpf\x)^(x) . (4.9) 
Note also that Eq. ([4.7| ) implies that the constant C can be written as 

/oo 

dxp ( °\x)v R (x) . (4.10) 
-oo 
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OC 



(ur) = - I dxxp^(x)w R (x) , (4.14) 



Thus the VEVs (z/j) (i = R, I) can be expressed solely in terms of the functions p\ (x) and 
<fii(x) {i = R, I), which can be calculated by standard Monte Carlo simulations. Clearly 
the derivation presented here is quite general, and we can obtain similar formulae in any 
complex action system. 

In the present case, we can slightly simplify the formulae due to the symmetry under 
( p.6|) . Using the properties 

M^T = <Pr(x) , (4-11) 
ipxixf = tfii-x) , (4.12) 
pf\-x) = p[°\x), (4.13) 

which follow from the symmetry, we arrive at 

,,<«), 

C = I dxp$\x)w R (x) , (4.16) 

J — OO 

where the weight functions Wi(x) are defined by 

w R {x) = (cosr) RiX , (4.17) 
wx(x) = (sinr) I)X . (4.18) 

Thus the problem reduces to the calculation of the four real functions pf^ (x), Wi(x) (i = R, I), 
which we will now discuss. 

4.2 Monte Carlo evaluation of pf\x) and Wi(x) 

In order to obtain Wi(x) (i = R, I), we need to simulate (|4.4|) . In practice, we simulate a 
partition function where the 5-function is replaced by a sharply peaked potential 



-oo 

2i ro ° 



(ui) = — / dxxpi' (x)wi(x) , (4-15) 
C Jo 



df e" So e- w . (4.19) 



In this study we use a Gaussian potential 

V(x) = ^(x-0 2 , (4.20) 

where 7 and £ are real parameters. The results are insensitive to the choice of 7 as far as 
they are big enough (we used 7 = 1000.0). Let us denote the VEV associated with the 
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partition function ( |4.19| ) as (0)iy. The expectation value (cosr)^ represents the value of 
wr(x) at x = (^r)r,v, while the expectation value (shiT)^ represents the value of wi(x) at 
x = (vi)i,v- 

In fact we can also obtain the functions pf\x) from the same simulation (|4.19|) . For that 
we note that the distribution of Vi in the system (}4.19|) is given by 



Pl y[x) d ^ (8(x - Vi ))iy oc P f\x) e~ v ^ 



(0), 



(4.21; 



which typically has a peak. The position of the peak, which we denote as x = x, is given by 
the solution to 







d 

dx 



\n Pl (x) = fi°\x)-V'(x) 



where we have introduced 



{0) <x)^ j-\npf\x) 



(4.22) 



(4.23) 



Therefore, the quantity V'(x) represents the value of f-°\x) at x = x. Since the value of 
7 is taken to be large, the peak is sharp and we can safely approximate x by the VEV 
(vi)i,v- Once we obtain f\ (x) for various x, we can calculate pf\x) by integrating ( |4.23| ) 
and exponentiating, where the integration constant can be determined by the normalization 
of the distribution pf\x). If the distribution pf\x) would be Gaussian (as it is the case 
near its peak; see Appendix), the above procedure is exact even for finite 7. 

Monte Carlo simulation of ( (4.19| ) can be performed by using the Hybrid Monte Carlo 
algorithm in much the same way as in Refs. fTl ). The required computational effort for the 
present model is 0(iV 3 ). 



4.3 The virtues of the method 

Comparing (|4.14| ), ( |4 . 1 5| ) and ( |4. 16| ) with Eq. p.8j ), we notice that 

/oo 
dxx p^\x) %(x) , 
-00 

POO 

(z/ I sinr) = 2 dx x pj°\x) wi(x) , 
Jo 

POO 

(cosr)o = / dx p^ (x) w R (x) . 



(4.24) 
(4.25) 
(4.26) 



Thus the new method as it stands simply amounts to using the standard reweighting formula 
([I8D, but calculating each VEV by using ( [4.24|) , (|4.2q) and (|4.26|) . We now explain the 
virtues of the present method. 

If we are to obtain the VEVs by directly simulating the system ( |3.5| ), for most of the 
time we sample configurations whose Uj takes a value close to the peak of pf\x). However, 
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from the r.h.s. of fl4.24| )~( |4.26| ) it is clear that we have to sample configurations whose U{ 
takes a value where pf\x)\wi(x)\ becomes large, in order to obtain the VEVs accurately. In 
general these two regions of configuration space have little overlap, which makes the sampling 
ineffective. This is the overlap problem. Since the overlap becomes exponentially small as the 
system size increases, this composes some portion of the complex action problem. The use 
of ( 4.24j) ~ (|4.26|) and calculating the relevant functions as explained in the previous section 
avoids this problem by 'forcing' the simulation to sample the important region. 

The knowledge of the weight factor Wi(x), which is provided by the present approach, 
allows us to probe directly the effect of the imaginary part F on the observable of our concern. 
We can understand which values of the observable are enhanced or suppressed by the effect of 
T. On the other hand, the standard reweighting technique simply gives the integrals on the 
l.h.s. of ( f4.24j) ~( |4.26|) , from which we can hardly imagine how they resulted from the effect 
of T. In the particular model we are studying, the imaginary part T has a dramatic effect on 
the VEV of the quark number as we discussed in Section [| We will see in the next section 
that the weight factor Wi(x) indeed provides a clear understanding of this phenomenon. 

As we also see in the next section, a Monte Carlo calculation of the weight factor Wi(x) 
becomes increasingly difficult with the system size. In that sense the complex action problem 
is still there. This should be contrasted with the meron-cluster algorithm |l5j , which has been 
applied to a special class of complex action systems with computational efforts increasing at 
most by some power of the system size. However, as mentioned, the factorization method 
eliminates the overlap problem and this is a substantial step forward which allows us to get 
closer to the thermodynamic limit. Indeed, we are able to obtain the thermodynamic limit 
of the 'quark number density' in this random matrix model with modest computer resources. 
See also Ref. |16]| for an idea to ameliorate the overlap problem in 'real' QCD at finite baryon 
density by interpolations in the (/i, T) plane (notice however [|17j1 ). 

Using the generic scaling properties of the weight factor Wt(x), one may extrapolate the 
results obtained by direct Monte Carlo evaluations to larger system size. Such an extrap- 
olation is expected to be particularly useful in cases where the distribution function turns 
out to be positive definite. In those cases we can actually even avoid using the reweighting 
formula ( |3.8| ) by reducing the question of obtaining the expectation value to that of finding 
the minimum of the free energy, which is (minus) the log of the distribution function. Here, 
the error in obtaining the scaling function propagates to the final result without significant 
magnifications. Therefore, the extrapolation can be a powerful tool to probe the thermody- 
namic limit from the accessible system size. Indeed such a technique has been used in Ref. 
U to discuss the dynamical generation of four- dimensional space-time in a nonperturbative 
formulation |18| of type IIB superstring theory in ten dimensions. 
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5 Reproducing exact results by the new method 



In this section we apply the factorization method to the model Q2.1|) as described in the 
previous section and show that it reproduces the known results (v) given by (|2.8|) for arbitrary 
N. Let us recall the dramatic difference between ( |2.9|) and (|3.9|) . It is particularly interesting 
to see how this occurs due to the effects of T in the present approach. First we focus on 
two values of /x, /x = 0.2 and /x = 1.0, which are on opposite sides of the first order phase 
transition point /x = /x c = 0.527 • ■ -. 





N 




i (vi) 


<"> 


(u) (exact) 


0.2 


8 


0.0056(6) 


-0.1970(5) 


-0.1915(7) 


-0.20000. . . 


0.2 


16 


0.0060(4) 


-0.1905(13) 


-0.1845(13) 


-0.20000. . . 


0.2 


24 


0.0076(9) 


-0.1972(14) 


-0.1896(17) 


-0.20000. . . 


0.2 


32 


0.0021(8) 


-0.1947(19) 


-0.1927(25) 


-0.20000. . . 


0.2 


48 


0.0086(37) 


-0.2086(54) 


-0.2000(88) 


-0.20000. . . 


1.0 


8 


0.8617(10) 


0.1981(13) 


1.0598(12) 


1.066501... 


1.0 


16 


0.8936(2) 


0.1353(6) 


1.0289(5) 


1.032240. . . 


1.0 


32 


0.9207(1) 


0.0945(2) 


1.0152(3) 


1.015871... 



Table 1: Results of the analysis of (u) described in the text. Statistical errors computed by 
the jackknife method are shown. The last column represents the exact result (|2.8j ) for (u) at 
each /x and N. For /x = 0.2 eq. ( |2.8| ) yields (u) = —0.2 with an accuracy better than 1 part 
in HT 9 . 



5.1 fi < fi c 

We start with the results for /x = 0.2. In Figs. || and [3], we plot and wi(x) respectively. 

In Figs. |] and [5], we show the results for the functions f^°\x). (The behavior of these 
functions can be understood theoretically as discussed in the Appendix.) By integrating 
these functions and exponentiating, we obtain the pf\x), which is shown in Figs. ^| and [7| 
In Fig. H and [5], we show and pf'\x)wi(x), respectively. Using them we obtain 

(z/r) and i {u\). Summing these values, we get (u), which should be compared with the exact 
result obtained from Eq. ( |2.8| ). The results are shown in Table [I]. 

Note that the sign change of xx>r(x) occurs near the peak of p^\x), so that the product 
Pr\x)wh(x) has a positive regime and a negative regime, which cancel each other resulting 
in (z/r) ~ 0. Thus the main contribution to (u) comes from the imaginary part (y{) in 
contrast to the results (|3.10| ) for the phase quenched system. This is consistent with the 
theoretical argument at the end of Section ||| that the contribution to (v) comes solely from 
the imaginary part in the large A" limit. 
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Figure 2: The weight factor w R (x) is plotted against x for iV = 8, 16, 24, 32, 48, 64 at jj, = 0.2. 
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Figure 3: The weight factor wi(x) is plotted against x for N = 8, 16, 24, 32, 48, 64 at \i = 0.2. 
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solid line represents the asymptotic behavior (|A.4j) discussed in the Appendix. 
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solid line represents the asymptotic behavior ( |A.5|) discussed in the Appendix. 
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position of the peak is approaching the large N result (z/r)o = 0.2 for the phase quenched 
system. 



6 







8 






16 






24 






32 






48 






64 






96 















0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

x 

Figure 7: The function pj 0) (a;) is plotted for N = 8, 16, 24, 32, 48, 64, 96 at p = 0.2. 
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5.2 fi > fi c 

We perform a similar analysis at p = 1.0. The results are presented in Figs. ITDL . . . , 17 



Again the results nicely reproduce the exact results as can be seen in Table [l]. Note in 
particular that the finite N effect is 1.6% at iV = 32, meaning that the accessible values of 
N are large enough to extract the large N limit. 

Here, wn(x) is approximately constant in the region where p^ (x) is peaked, so the 
shape of the product p^ (x)wn(x) is similar to p& (x). On the other hand, the peak of 
Pi(x) at x = is slightly shifted by multiplying wj(x), but the first moment of the product 
Pi \x)wi(x) is still small. This is also the case for p = 0.2, but the difference comes from 
the normalization constant C in the formula ( |4.15|) for (y\). As in (|4.16|) , the constant 
C is obtained by integrating p^(i)wr(i), where cancellations occur at p = 0.2, but not 
at p — 1.0. As a result we obtain (b>i) ~ at p, — 1.0 as N increases. Thus the main 
contribution to (u) comes from the real part (^r), and moreover, it is close to (^r)o- Again 
this is consistent with the theoretical argument given at the end of Section || 

It is interesting that the changes from positive to negative for p = 0.2, but it 

changes from negative to positive for p = 1.0. Similarly w\{x) is positive at x > for 
p = 0.2, but it is negative at x > for p = 1.0. Thus the behavior of Wi(x) changes 
drastically as the chemical potential p crosses its critical value p c . 

5.3 fl ~ fl c 

Let us see in more detail what is going on in the critical regime. In Figs. |TB| and |TJ] we 
plot and Wi(x) respectively for iV = 8 at various p. The final results for (z/j) are 

plotted in Fig. ^0] and the corresponding data are listed in Table |[ Note that we were able 
to reproduce the exact result even at the turning point (p = 0.6139). These results provide 
a clear understanding of how the first order phase transition occurs due to the effects of T. 

The fact that \wi(x) | becomes small near the critical regime reveals an increasing difficulty 
in approaching the critical point. In Figs. |2l], |22|, we plot ln(max x |u>i(a;)|) against \p — p c \ 
at N = 16, 32 for p < p c and p > p c respectively. Our data can be nicely fitted to 

ln(max = — aexp(— b\p — p c \) + const. . (5.1) 

X 

For fixed p, on the other hand, we observe that max x |wi(x)| decreases exponentially as 



oc exp(— cN). We plot the result for p = 0.2 in Fig. [23|. This is consistent with the scaling 



behavior of the weight function, which plays a crucial role in the extrapolation [|TJ as we 



mentioned at the end of Section 4.3 
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Figure 10: The weight factor w R (x) is plotted against x for N = 8, 16, 32 at ji — 1.0. 
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Figure 11: The weight factor tui(:r) is plotted against x for iV = 8, 16, 32 at /i = 1.0. 
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Figure 12: The function f^\x) is plotted for iV = 8,16,32 at fi — 1.0. The solid line 
represents the asymptotic behavior ( |A.4|) discussed in the Appendix. 



5 




-5 
-10 
-15 
-20 
-25 



m. X ^ 5( W 




V 



8 
16 
32 
-4/x 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

x 

Figure 13: The function fi(x) is plotted for iV = 8,16,32 at /i — 1.0. The solid line 
represents the asymptotic behavior ( |A.5| ) discussed in the Appendix. 
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Figure 14: The function p^\x) is plotted for N = 8, 16, 32 at n = 1.0. The position of the 
peak is approaching the large N result (z/r)o = 1 for the phase quenched system. 
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Figure 18: The weight factor is plotted against x for N = 8 at various \x. The behavior 

changes drastically as \i crosses the critical point. 
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Figure 19: The weight factor w\{x) is plotted against x for N = 8 at various /j,. The behavior 
changes drastically as \i crosses the critical point. 
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Figure 20: The VEV (v) is obtained by the factorization method at iV = 8 for various \i 
including the critical regime. Statistical errors computed by the jackknife method are also 
shown. The dashed line represents the exact result fl2.8|) for (v) at N = 8. 



X 





-0.02 
-0.04 
-0.06 
-0.08 
-0.1 
-0.12 
-0.14 
-0.16 
-0.18 




0.1 



0.15 



0.2 



0.25 

IH-Mcrl 



0.3 



0.35 



0.4 



Figure 21: The result of lnfmax-j. |u>i(x)|) is plotted against \fx — fj, c \ ior N = 16,32 at ji < fi c 
The lines are the fits to the behavior ( |5.1[ ). 
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The lines are the fits to the behavior ( |5.1|) . 
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Figure 23: The result of ln(max x \wi(x) |) at [A — 0.2 is plotted against N. A linearly 
decreasing behavior is observed at large N. 
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0.7 
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0.480(4) 


1.451(6) 


1.45516 
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0.444(3) 


1.417(4) 


1.42899 


0.8 
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0.949(1) 


0.390(2) 


1.339(3) 


1.34792 
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0.9054(3) 


0.2779(4) 


1.1833(6) 


1.19253 


1.0 


8 


0.8622(1) 


0.1999(2) 


1.0620(2) 


1.06650 



Table 2: Results of the analysis of (y) for N = 8 at various fi including the critical regime. 
Statistical errors computed by the jackknife method are also shown. The last column repre- 
sent the exact result ( |2.8|) for (u) at iV = 8. 

6 Summary 

In this article we have clarified the properties of the factorization method for systems with 
a complex action. This method circumvents the overlap problem, and we hope that this 
feature alone will make various interesting questions accessible by the present day computer 
resources. Indeed, we were able to reproduce the exact results for the quark number density 
in a schematic model for QCD at finite baryon density. The achieved system size was already 
large enough to obtain the thermodynamic limit. We therefore expect that the factorization 
method is useful to explore the phase diagram of the 'real' finite density QCD. The method 
itself is quite general, and it can be applied to any system with a complex action, although 
the actual gains may depend on the system. 

We also emphasize that in the case where the distribution functions turn out to be 
positive definite, the method can be even more powerful by utilizing the scaling property 
of the weight factor. This extrapolation appeared particularly valuable in the study of 
spontaneous symmetry breaking in the type IIB matrix model. Although this is a very 
interesting problem which would tell us a great deal about nonperturbative string dynamics 
and the dynamical origin of the space-time dimensionality, we expect that this method will 
be useful for other systems as well. 
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Appendix: Large N behavior of pf\x 

In this Appendix, we discuss the large N behavior of functions pf\x). From Figures 
[|, ^, O, [IB], we find that p[ (x) is well approximated by the Gaussian distribution near 
the peak, but there is a transition to a power-like tail (pf\x) oc x~ 4 ) at large \x\. The 
function fi(x) scales in this power regime. In the Gaussian regime, on the other hand, the 



function Tffi°\x), with the normalization factor 1/N, scales as is shown in Figs. |24| and |25 



for fi = 0.2. As a result the extent of the Gaussian regime shrinks as 1/yN. Such a behavior 
is obtained if N independently distributed eigenvalues of W contribute to (v). Indeed, the 
correlations of the eigenvalues of the matrix W in our model decay exponentially on the 
scale of the average level spacing. 

The observed power tail can be understood as follows. The point is that large values of 
the baryon number result from eigenvalues of the matrix W close to ±i/i. The probability of 
finding one eigenvalue A inside a radius a can be easily obtained for /i = from generalizing a 
calculation in the book of Mehta [19|] (the chapter on RMT's without hermiticity conditions). 
One finds that 

N 2 A 

P(3A < a) = —a 4 for a -> . (A.l) 
If we assume stationarity of this result we obtain 



P(\u\>x) = P[B\,\\-iu\< = —r. (A.2) 

1 '' Pl 2Nx 32N 2 x 4 V ; 



For the probability density we thus find 



f(M = *> = ass? • (A ' 3) 
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Figure 24: The function ±f£ ] (x) is plotted for \x = 0.2. A clear scaling behavior is seen in 
the linear regime, where the function crosses zero. 
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Figure 25: The function jjfi°\x) is plotted for \x = 0.2. A clear scaling behavior is seen in 
the linear regime, where the function crosses zero. 
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so that the real and imaginary parts of the quark number density behave as l/x A for large 
\x\. Since the level spacing distribution is a universal feature of Random Matrix Theories, 
we expect that such power-like tails are generic and also occur in QCD. 

In fact our results in Figures limn suggest that f^\x) and f$°\x) are given at 
large \x\ by 

fL°\x) ~ --±- (A.4) 

fi%) ~ --, (A.5) 
x 

which agrees with the above argument. 
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